module H2MolInfo
use prec
use iounits
use waterinfo
use regioninfo
implicit none
   !----------------------------------------------------------------------------
   private
   public  H2Mol_Init, Add_H2Mol, LH2Mol, H2_Init
   public  nH2atom, nH2Bond, H2Mass, H2Name
   public  H2MolID, H2Moltype, H2MolBond, H2MolPos, H2MolChg
   public  HHBondCoeff, HHPair
   !----------------------------------------------------------------------------
   logical :: LH2Mol
   integer :: nH2Mol = 0, nH2atom = 0, nH2Bond = 0, nH2Angle = 0
   !----------------------------------------------------------------------------
   integer, allocatable :: H2MolType(:), H2MolBond(:,:)
   integer, allocatable :: H2MolID(:)
   real(q), allocatable :: H2MolPos(:,:), H2MolChg(:)
   !----------------------------------------------------------------------------
   integer :: H2Type(2), H2Bond(2,1)
   real(q) :: H2Pos(3,2), H2Chg(2)
   real(q) :: H2Mass(1) = (/ 1.00794D0 /)
   !----------------------------------------------------------------------------
   ! J. Chem Phys 125, 014702, 2006; J Chem Phys 121, 10096, 2004. TIP3P-Ew
   real(q), parameter :: HHBondCoeff(2) = (/ 390.26D0, 0.7375D0  /)
   real(q), parameter :: HHPair(2,1)   = reshape((/0.2521D0,2.645D0/), (/2,1/))
   character(len=2)   :: H2Name(1) = (/ "H" /)
   !----------------------------------------------------------------------------
contains
   !
   subroutine H2Mol_Init()
   implicit none
      !-------------------------------------------------------------------------
      character(len=256) :: oneline
      integer            :: ioerr
      !-------------------------------------------------------------------------
      write(stdout, '(/,10x,"Would you like to insert H2 molecules into the system? (y/n)[n]: ", $)')
      read(stdin,'(A)',iostat=ioerr) oneline
      LH2Mol = ioerr.eq.0.and.(oneline.eq.'y'.or.oneline.eq.'Y')
      !
      if (.not.LH2Mol) return
      !
      ! initialize H2 molecule info
      H2Type      = (/ 1, 1 /)
      H2Bond(:,1) = (/ 1, 2 /)
      H2Chg       = (/ 0.D0, 0.D0 /)
      H2Pos(:,1)  = (/ 0.D0, 0.D0, 0.D0 /)
      H2Pos(:,2)  = (/ 0.74D0, 0.D0, 0.D0 /)
      !
      ! ask for total number of H2 molecules
      write(stdout, '(10x,"How many H2 molecules would you like to insert? : ",$)')
      nH2Mol = 0
      do while (ioerr.ne.0.or.nH2Mol.lt.1)
         read(stdin,*,iostat=ioerr) nH2Mol
      enddo
      nH2Atom = nH2Mol * 2
      nH2Bond = nH2Mol
      nH2Angle= 0
      !
      if ( allocated( H2MolType ) ) deallocate( H2MolType )
      if ( allocated( H2MolBond ) ) deallocate( H2MolBond )
      if ( allocated( H2MolPos  ) ) deallocate( H2MolPos  )
      if ( allocated( H2MolChg  ) ) deallocate( H2MolChg  )
      if ( allocated( H2MolID   ) ) deallocate( H2MolID   )
      allocate( H2MolType(nH2Atom), H2MolBond(2,nH2Bond), H2MolPos(3,nH2Atom), H2MolChg(nH2Atom), H2MolID(nH2Atom) )
      !
      call init_region()
      !
   return
   end subroutine

   subroutine Add_H2Mol()
   use cell, only: box, zboxlow, zboxhigh
   use wallinfo, only: LAddWalls
   implicit none
      !-------------------------------------------------------------------------
      integer :: i, istr, iend, iprev
      real(q) :: rndnum(6), TPI, rotatedH2(3,2)
      logical :: loverlap
      !-------------------------------------------------------------------------
      TPI  = 8.D0 * atan(1.D0)
      !
      do i = 1, nH2Mol
         loverlap = .true.
         istr     = (i-1) * 2 + 1
         iend     = i     * 2
         iprev    = istr  - 1
         ! To rotate and translate the molecule randomly, make sure that no atom will overlap.
         do while ( loverlap )
            call RANDOM_NUMBER( rndnum )
            rndnum(1:3) = rndnum(1:3) * TPI
            call rotate_molecular( rndnum(1), rndnum(2), rndnum(3), H2Pos, rotatedH2, 2)
            call translate_molecular( rndnum(4), rndnum(5), rndnum(6), rotatedH2, 2)
            if ( lregion ) loverlap = check_out_region(rotatedH2)
            if (.not.loverlap.and.LAddWalls ) then
               loverlap = any(rotatedH2(3,:).lt.zboxlow).or.any(rotatedH2(3,:).gt.zboxhigh)
            endif
            if (.not.loverlap) loverlap = check_H2Mol_position(iprev, rotatedH2)
         enddo
         ! position, atom type and chage of the i-th molecule
         H2MolPos(:, istr:iend)  = rotatedH2
         H2MolType(istr:iend)    = H2type
         H2MolChg(istr:iend)     = H2Chg
         H2MolID(istr:iend)      = i
         !
         ! bonds of the i-th molecule
         H2MolBond(:,i) = H2Bond(:,1) + iprev
         !
      enddo
   return
   end subroutine

   logical function check_H2Mol_position(iprev, pos)
   use cell, only: atpos, natom, box, rbox
   use wallinfo, only: LAddWalls, WallPos, nWallAtom, WallRemain
   implicit none
      !-------------------------------------------------------------------------
      integer, intent(in) :: iprev
      real(q), intent(in) :: pos(3,2)
      !-------------------------------------------------------------------------
      real(q) :: Rmin2(2) = (/ 1.0D0, 1.0D0 /), Rwmin = 2.0D0
      !-------------------------------------------------------------------------
      integer :: i, j, k
      real(q) :: R(3), dist
      !-------------------------------------------------------------------------
      check_H2Mol_position = .false.
      !
      ! all liquid atoms
liquid:do i = 1, natom
         do j = 1, 2
            R = pos(:,j) - atpos(:,i)
            R = R - nint(R*rbox)*box
            dist = R(1)*R(1) + R(2)*R(2) + R(3)*R(3)
            check_H2Mol_position = dist.lt.Rmin2(j)
            if (check_H2Mol_position) exit liquid
         enddo
      enddo liquid
      ! Wall atoms
      if ( LAddWalls.and.(.not.check_H2Mol_position) ) then
wall:    do i = 1, nWallAtom
            if ( WallRemain(i).eq.0 ) cycle
            do j = 1, 2
               R = pos(:,j) - WallPos(:,i)
               R(1:2) = R(1:2) - nint( R(1:2)*rbox(1:2) ) * box(1:2)
               dist   = R(1)*R(1) + R(2)*R(2) + R(3)*R(3)
               !
               check_H2Mol_position = dist.lt.Rwmin
               if (check_H2Mol_position) exit wall
            enddo
         enddo wall
      endif
      ! Water atoms
      if ( lwater.and.(.not.check_H2Mol_position) ) then
water:   do i = 1, nH2Oatom
            do j = 1, 2
               R = pos(:,j) - waterpos(:, i)
               R = R - nint( R*rbox ) * box
               dist = R(1)*R(1) + R(2)*R(2) + R(3)*R(3)
               check_H2Mol_position = dist.lt.Rmin2(j)
               if (check_H2Mol_position) exit water
            enddo
         enddo water
      endif
      ! Previous H2 atoms
      if (.not.check_H2Mol_position) then
h2mol:   do i = 1, iprev
            do j = 1, 2
               R = pos(:,j) - H2MolPos(:, i)
               R = R - nint( R*rbox ) * box
               dist = R(1)*R(1) + R(2)*R(2) + R(3)*R(3)
               check_H2Mol_position = dist.lt.Rmin2(j)
               if (check_H2Mol_position) exit h2mol
            enddo
         enddo h2mol
      endif
   return
   end function check_H2Mol_position

   subroutine H2_Init()
   use cell
   use general
   implicit none
      !----------------------------------------------------------------------------
      integer  :: i, moltype
      !----------------------------------------------------------------------------
      subnamep = subname
      subname  = 'H2_init'
      !
      if ( allocated( bondsMol  ) ) deallocate( bondsMol  )
      if ( allocated( anglsMol  ) ) deallocate( anglsMol  )
      if ( allocated( dihrsMol  ) ) deallocate( dihrsMol  )
      if ( allocated( attypMol  ) ) deallocate( attypMol  )
      if ( allocated( atposMol  ) ) deallocate( atposMol  )
      if ( allocated( molIDMol  ) ) deallocate( molIDMol  )
      if ( allocated( rotatedMol) ) deallocate( rotatedMol)
      if ( allocated( chargeMol ) ) deallocate( chargeMol )
      if ( allocated( EMass     ) ) deallocate( EMass     )
      if ( allocated( sigma     ) ) deallocate( sigma     )
      if ( allocated( epison    ) ) deallocate( epison    )
      if ( allocated( bondcoef  ) ) deallocate( bondcoef  )
      if ( allocated( anglcoef  ) ) deallocate( anglcoef  )
      if ( allocated( torscoef  ) ) deallocate( torscoef  )
      if ( allocated( DistMin   ) ) deallocate( DistMin   )
      if ( allocated( EName     ) ) deallocate( EName     )
      !
      ntypeMol   = 1  ! 1: H
      nbndtypMol = 1
      nangtypMol = 0
      ndihtypMol = 0
      !
      allocate( EMass(ntypeMol), sigma(ntypeMol), epison(ntypeMol),           &
      &         bondcoef(2, nbndtypMol), DistMin(ntypeMol), EName(ntypeMol+1) )
      EMass    = (/ 1.00794D0 /)
      sigma    = (/ 2.645D0   /)
      epison   = (/ 0.2521D0  /)
      DistMin  = (/ 1.0D0     /)
      EName    = (/ 'H'       /)
      !
      bondcoef = reshape( (/ 390.26D0, 0.7375D0 /),(/2, nbndtypMol/) ) ! H-H
      !
      molform    = 'H2'
      molname    = 'H2'
      den0       = 1.D0
      natomMol   = 2
      nbondsMol  = 1
      nanglsMol  = 0
      ndihrsMol  = 0
      !
      allocate( bondsMol(3, nbondsMol), attypMol(natomMol),                   &
      &         atposMol(3, natomMol), chargeMol(natomMol), molIDMol(natomMol) )                
      !
      attypMol = (/ 1, 1 /)
      !
      atposMol(:, 1) = (/ 0.D0, 0.D0, 0.D0 /)
      atposMol(:, 2) = (/ 0.74D0, 0.D0, 0.D0 /)
      chargeMol      = (/ 0.D0, 0.D0 /)
      !                    Type, atom 1, atom 2
      bondsMol(:, 1) = (/   1,  1,  2 /)
      molIDMol       = (/ 1, 1 /)
      !
      nMol      = 1
      molmass   = 0.D0
      do i = 1, natomMol
         molmass = molmass + EMass( attypMol(i) )
      enddo
      forall ( i=1:3 ) minbox(i) = maxval(atposMol(i,:)) - minval(atposMol(i,:))
      minlen = maxval( minbox )
      allocate( rotatedMol( 3, natomMol) )
      !
      subname = subnamep
   return
   end subroutine
end module
